Simulation Details
Suppose
N vibrating molecules (0, 1, 2,...,N-1) modelled as harmonic
oscillators. Their quantized energy levels, i, are given by:
εi
= (i + 1/2) hν; (i = 0, 1, 2,...)
where
h is Planck's constant and ν the oscillator frequency.
The
energy levels are equally spaced (εi+1- εi
= hν) and taking the zero-point energy (1/2 hν) as a reference, there is
no ambiguity in denoting the level energies by 0, 1, 2,... That is,
energy "0" corresponds to 1/2 hν, energy "1" to hν, energy "2" to 2 hν,
and so forth.
In
this scheme (see Variables panel):
Integer
array npi [ ] contains the number of molecules at the energy levels :
npi[0]
- number of molecules at level 0; npi[1] - number of molecules at level 1;
...
Integer
array en [ ] contains the energies of molecules i:
en[0]
- energy of molecule 0; en[1] - energy of molecule 1; ...
Note:
as the energies of any two successive levels differ by 1u.e. (hν), the
energy of a molecule at level i, is equal to that level label, i. For
example, if molecule 30 is at level 20, then en[30] = 20. Next, it is seen
that proceeding in integer arithmetic, for manipulating those arrays,
facilitates the programming.
a)
Initially, 1000 molecules are set at level 10 (see Initialization
panel). Therefore, the total energy of the system is 1000*10.
b)
Two molecules, i1 and i2, are chosen at random (see method collide( ) in
the Custom panel). If the chosen molecules are different, then they
collide: one molecule increases the energy by 1 u.e., the other decreases
the energy by the same amount, and the array storing the number of
molecules at each energy level is updated accordingly. Therefore, the
total energy is kept constant. Note, as suggested above, how the code
takes advantage of the integer arrays and that the energy array (en), is
the index of the number of particles array (npi).
For
example, if the molecule 30 at level 20 (en[30] = 20) jumps to level 21,
then:
npi[
en[30] ] = npi [ en[30] ] -1 (the number of molecules at level 20
decreases 1 unit);
en[30]
= en[30] +1 (en[30] = 21, i.e, the molecule is now at level 21);
npi[
en[30] ] = npi[ en[30] ] + 1 (level 21 is updated);
Calculation
of ln W (see methods entropy( ) and factorial ( int K) in the
Custom panel)
The
basic expressions are:
ln
W = ln N! - ∑ ln ni!
ln
N! = ln 1 + ln 2 +...+ ln N; similar for the other ln's
Calculation
of β for the theoretical Boltzmann's distribution (see
Initialization panel)
This
calculation is a bit tricky. Let´s define:
x
= exp (-β) where 0<x<1
For
the present particular case the average energy per molecule is:
1/N
∑ni εi = 1/Z ∑εi x εi
= 10 u.e,
since,
initially, the N molecules were set at level 10, and the total energy is
kept constant. Z is the partition function ( ∑xεi ).
For
the present initial conditions, preliminary tests have shown that the
energy levels worth to consider are the first 40 ones. Thus, expanding the
last equation to the first 41 levels:
x
(9 + 8x + 7x2 + ... + x8 - x10 - 2x11
- ... - 31x40) = 0
The
numerical root (between 0 and 1) of the algebraic equation of degree 40,
is:
x =
0.907 whence β = 0.10, the value used to work out the theoretical
Boltzmann's distribution.
In
order to prevent out-of-bound errors, a maximum number of 200 harmonic
oscillator levels have been considered (though only the first 41 are used
in the calculations and graphical representation) since, every now and
then, molecules can jump to higher levels.
However,
if other initial conditions are chosen, care must be taken of the
appropriate number of energy levels, the recalculation of β and the arrays
dimension.